Packages you will need:

library(RColorBrewer)
library(leaflet)
library(ggplot2)
library(mapview)
library(raster)
library(htmltools)
library(sf)

Read in your data

There are many ways to read spatial data into R
Here we get admin level 2 boundaries for Nepal from GADM (using raster package) and then convert the objcect to an sf class.
If you don’t have internet, you can read in shapefiles from your drive (as you might do in ArcGIS). R-packages also come with inbuilt spatial datasets- see this Mapview tutorial which uses an inbuilt dataset on breweries in a region of Germany

Nepal2 <- getData('GADM', country='NPL', level=2)

Nepal2 <- st_as_sf(Nepal2)  #  convert to an sf object  (for later processing)

1. Basic, fast plotting with Mapview

mapview(Nepal2)
# colour districts by name
mapview(Nepal2, zcol = "NAME_2")

Let’s generate some dummy data:

## generate random number for cases
Nepal2$cases <- sample(0:100, nrow(Nepal2), replace=T)

mapview(Nepal2, zcol = "cases")

It’s easy to define breaks using seq() function

mapview(Nepal2, zcol = "cases", at = seq(from = 1, to = 100, by =10), legend = TRUE)

Create a new field to show (made-up) rate

Nepal2$population <- sample(10000:30000, nrow(Nepal2), replace=T)
Nepal2$rate <- round(Nepal2$cases/Nepal2$population*10000,1)
mapview(Nepal2, zcol = "rate")
# generate some made up health facilities
facilities <- st_sample(x = Nepal2, size = 50, type = "random", exact = TRUE)


## plot Nepal districts and health facilities
mapview(Nepal2, zcol = "cases") + facilities
## change colour of facilities
mapview(Nepal2, zcol = "cases") + mapview(facilities, col.regions = "red")
# define your own palette to colour the areas 
pal <- colorRampPalette(c("white", "brown1"))

mapview(Nepal2, zcol = "cases", at = seq(0, 100, 20),  col.regions = pal(5), legend = TRUE)

2. More advanced mapping with Leaflet

Leaflet has many more options to customise the presentation of your map
The learning curve is steeper than for mapview (but there’s lots of online examples)
It’s ideal if you want to produce maps to present on a website

# Use leaflet to display a blank map with OpenStreetMap imagery 
map <-leaflet() %>%                           ## using the pipe operator to chain operations together: here the blank map is chained to 
  addProviderTiles(providers$OpenStreetMap)   ## the first leaflet option we are applying
map

Add Nepal districts to the OpenStreet base map

Map_Nepal<-map %>%                               
  addPolygons(data=Nepal2)
Map_Nepal

Define colour schemes for repeatable plotting

colorBin() is a leaflet function that maps continuous numerica data onto a palette First make a gradient with as many colours as you like e.g. red-yellow-green
Then produce a palette mapping the values in your dataframe to this gradient
Or you can map the data to one of R’s pre-defined palettes (e.g. “Blues”)

My_gradient <- c("white", "brown1" )       ## brown1 is actually red
palTotal <-  colorBin(My_gradient, domain = Nepal2$cases, bins = quantile(Nepal2$cases))  

palRate <-  colorBin("Blues", domain = Nepal2$rate, bins = quantile(Nepal2$rate))

# we will call these palettes in options when we make the map

Static map of number of cases by district

This generates a simple, non-interactive map
Here we use different options to define the outline and fill style, and assign the layer to a group
There are numerous other options for you to explore here For example, try making the lines dashed and white instead of solid and black

Map_Nepal <- map %>%                             ## use the pipe operator to add elements 
                                                 ## to the OpenStreet base map      
  addPolygons(data = Nepal2,                       ## Use data from the spatial polygons dataframe Nepal2
              color = 'black',                   ## line color    
              weight =2,                         ## line width
              opacity = 1,                       ## line opacity
              fillColor =  ~palTotal(cases),      ## fill polygons using the palette we generated
              fillOpacity = 1,                   ## opaque fill
              group = "Number of cases")         ## allows this element to be linked e.g. to legend           
Map_Nepal

Add layers control panel and a legend

When we start adding more layers to the map, we will need a control panel so the user can choose which layers to display
We will also need a layers panel that should be linked to the data and the control panel

Map_Nepal <- Map_Nepal %>%                       ## again we use the pipe operator 
                                                 ## to add complexity to our existing map object
  
## E.g. add the "Number of cases" group to the layers panel.   
  addLayersControl(overlayGroups = c("Number of cases")) %>%   
  
## Currently this group only contains  the polygons showing number of cases.
  
## Add a Legend  
  addLegend(pal = palTotal,                           ##  use same palette as for polygons     
            values = Nepal2$cases,                    ##  use same values 
            opacity = 1,                              ##  self-explanatory
            title = "Number of cases",                ##  "
            position = "bottomleft",                  ##  "
            group = "Number of cases",                ## links this legend with the polygons
            labFormat = labelFormat(digits = 0))      ## define decimal places     
Map_Nepal

Adding more layers

We can add layers showing other fields just by changing the fillColor argument and assigning it to a new group
To make your maps look consistent, its usually best to keep the other options (e.g. line style) the same

Map_Nepal <- Map_Nepal %>%
  addPolygons(data = Nepal2, 
              weight = 2, 
              opacity = 1,
              color = 'black',
              fillColor =  ~palRate(rate),     ## colour by rate
              fillOpacity = 1,
              group = "Rate") %>%             ## new group
  
# Layers control
addLayersControl(overlayGroups = c("Number of cases", "Rate"),            ## as before but with 2 groups
                   options = layersControlOptions(collapsed = FALSE))%>%  ## this keeps the panel expanded
# Legend
addLegend(pal = palRate, 
            values = Nepal2$rate, 
            opacity = 1, 
            title = "Rate (per 10,000)",
            position = "bottomleft", 
            group = "Rate",
            labFormat = labelFormat(digits = 1))

Map_Nepal

Add interaction: highlight and label districts on hover

This is a long chunk of code but it is just putting together what we have done above, and adding highlight and label options.
You could also look at this Example from leaflet

## define labels using the sprintf() function to format text
number_labels <- sprintf(
  "<strong>%s</strong><br/>%d cases",    # this is C-style formatting which means 
  Nepal2$NAME_2, Nepal2$cases            # "print NAME_2 in bold font, return, 
                                         #  print the number of cases in normal font"
) %>% lapply(htmltools::HTML)

rate_labels <- sprintf(
  "<strong>%s</strong><br/>%g cases per 10,000",
  Nepal2$NAME_2, round(Nepal2$rate, digits = 1)
) %>% lapply(htmltools::HTML)

## Add these labels to the map, along with highlighting options
## Number layer - inital options as above
Map_Nepal <- map %>%    
  addPolygons(data=Nepal2,   
              fillColor=  ~palTotal(cases),             
              fillOpacity = 1,                   
              color = 'black',                 
              weight =2,                         
              opacity = 1,                       
              group = "Number of cases",
              highlight = highlightOptions(             # When user hovers, give the polygon:
                weight = 4,                             # A thicker outline
                color = "#666"),                        # A grey outline
              label = number_labels,                    # The label we defined as number_labels    
              labelOptions = labelOptions(
                textsize = "15px")) %>%
  
  ## rate
  addPolygons(data=Nepal2, 
              fillColor=  ~palRate(rate),
              fillOpacity= 1,
              color= 'black',
              weight=2, 
              opacity= 1,
              group = "Rate",
              highlight = highlightOptions(
                weight = 4,
                color = "#666"),
              label = rate_labels,
              labelOptions = labelOptions(
                textsize = "15px")) %>%
  
  
  # Layers control
  addLayersControl(overlayGroups = c("Number of cases", "Rate"),  
                   options = layersControlOptions(collapsed = FALSE))%>% 
  
  # Legend
  addLegend(pal = palTotal,                   
            values = Nepal2$cases,                   
            opacity = 1,                              
            title = "Number of cases",               
            position = "bottomleft",                 
            group = "Number of cases",
            labFormat = labelFormat(digits = 0))%>% 
  
  addLegend(pal = palRate, 
            values = Nepal2$rate,
            opacity = 0.7, 
            title = "Rate (per 10,000)",
            position = "bottomleft", 
            group = "Rate",
            labFormat = labelFormat(digits = 1))%>% 

hideGroup("Rate")   ## This means only the first group will be displayed initially

Map_Nepal

Finally we customise other elements on the map

If you are making the map for a funder, a programme or any other external partner, you may need to add their logo
The url option is really useful if your maps are going to appear on a separate webpage

Map_Nepal <- Map_Nepal %>%  
                  addLogo(img = "http://blogs.lshtm.ac.uk/inthenews/files/2017/12/LSHTM-logo-bw-1.jpg",
                  url = "https://www.lshtm.ac.uk/",
                  position = "topleft",
                  width = 120,
                  height = 60) %>% addMeasure() 
                                     

Map_Nepal